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Spherical Gravitational Wave Detectors: 
an approach for real time data analysis 
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We have developed a full model to simulate spherical detectors where all main sources of noise 
are considered. We have built a computer code for determining the source direction and the wave 
polarization (solution of the inverse problem) in real time acquisition. The digital filter used is a 
simple bandpass filter. The "data" used for testing our code was simulated, which considers both 
the source signal and detector noise. The detector noise includes the antenna thermal, back action, 
phase noise, series noise and thermal from transducer coupled masses. The simulated noise takes 
into account all these noises and the antenna-transducers coupUng. The detector transfer function 
was calculated for a spherical antenna with six two-mode parametric transducers. From the results 
we determined that spherical detectors are able to locate an astrophysical source in the sky in real 
time as long as the signal-to- noise ratio for a burst signal is equal to or higher than 3. 
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I. INTRODUCTION 

Spherical gravitational wave detectors seem like an ex- 
cellent low cost alternative for GW detection. They have 
omnidirectional sensitivity and can give us some impor- 
tant information about an incident gravitational wave: 
its direction and tensorial components (polarization). 

Moreover a sphere can be thought as 5 detectors in 
a single instrument. The cross-correlation between nor- 
mal modes reduces noise effects and increases detection 
possibility. 

Two of these instruments are being developed: the 
Brazilian Mario SCHENBERG and the Dutch Mini- 
GRAIL Q and a third one is planned: the Sfera antenna 
0. Both detectors have similar features and their data 
analysis approach can be used in the same way. 

Robert Forward was the first to suggest a sphere as 
the antenna element in the early 1970's ^3]. Ashby and 
Dreitlein studied the reception of GWs by an elastic self- 
gravitating spherical antenna and found a set of equa- 
tions to treat such problem 4]. Following up this idea 
Wagoner and Paik found the lowest eigenvalues for the 
monopole and quadrupole modes of a uniform elastic 
sphere and calculated its cross-section 5j. 

In the 1990's, Johnson and Merkowitz studied the 
antenna-transducer coupling problem and found an op- 
timum configuration, which minimizes the number of 
transducers while keeping them in a symmetric distri- 
bution on the antenna: the truncated icosahedron (TI) 
configuration 6, 7]. Magalhaes and collaborators showed 
that distributions with more resonators can also have 
interesting symmetric properties 0, . Lobo also sug- 



gested a more elaborated symmetric distribution ^lOj . 
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FIG. 1: Schematic diagram of the Brazilian spherical GW 
detector Schenberg. 

For more than a decade financial and technologi- 
cal difficulties had stalled the construction of spherical 
GW detectors but nowadays some research groups are 
working on this. One of these is the Brazilian Mario 
SCHENBERG detector (see Figure [ij. The SCHEN- 
BERG has a 65 cm-diameter and ~ 1150 kg copper alloy 
[Cu(94%)Al(6%)] spherical antenna which is suspended 
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by vibrational isolators. It will be kept under tempera- 
tures below O.IK inside cryogenics chambers 11]. These 
features are shared with the Dutch MiniGRAIL detector. 
The main difference between both is that the Brazilian 
detector will use parametric transducers while the Dutch 
data will be taken from capacitive transducers. 

MiniGRAIL is ahead in the development. It has al- 
ready started commissioned operations and soon will 
start scientific runs. Even though SCHENBERG is only 
now ready to start a commissioned phase at 4.2 K, it has 
a data acquisition system and a data analysis software al- 
most ready for taking scientific data. We have developed 
a model based on harmonic oscillators to simulate their 
behavior. The development of software, which outputs 
simulated data, allowed us to create a methodology for 
initial data treatment. These procedures are presented 
in this paper along with how we can obtain information 
about GW sources directly from the detector data. 



II. THE MODEL AND THE SIMULATED DATA 

In order to develop a model that represents the re- 
sponse of a spherical detector to gravitational wave sig- 
nals we created a set of displacement equations for the 
antenna plus resonator system. These equations (based 
on the linear elastic theory) are ruled by physical param- 
eters such as element masses, spring constants, dumping 
coefficients and all coupling effects from electronic de- 
vices. The model considers the dynamics of a spherical 
antenna coupled to a set of electro-mechanical transduc- 
ers and supply voltage pre- amplifiers readouts. In other 
words the model output is equivalent to the data which 
will be taken from the detector (voltage measurements 
which generate a time series). 

The system dynamics can be expressed as a group of 
equations which represents coupled harmonic oscillators 
and is represented (already normalized) by 

w(t) + Hw(t) + Kw(t) = M-^WF(t), (1) 

where M, H and K are respectively: the mass, dumping 
and spring matrixes. w(t) is the displacement vector and 
F{t) is the forces vector. W denotes the forces transfered 
to the system. This kind of equation is easily solved when 
one uses its orthonormal form in frequency domain [T^ . 
The solution is the following 

w(cj) = UJ(cj)U'^M-^WF(cj), (2) 

where U is the eigenvectors matrix and U''' is its hermi- 
tian conjugated. The system response function is rep- 
resented by the matrix J~-^(cc;) = — cj^I + iujU^HXJ + D 
which couples all modes. Resonant frequencies are ruled 
by the diagonal matrix D (the eigenvalues of the system). 
The size of this set of equations depends on the number 
of transducers and its modes. 



A. The TI Configuration and the Mode Channels 

In order to increase the bandwidth we opted to use 
six 2-mode mechanical resonators coupled to the an- 
tenna surface according to the arrangement suggested by 
Johnson and Merkowitz: the truncated icosahedron (TI) 
configuration j^, . This symmetric distribution places 
transducers in positions which correspond to the centers 
of pentagonal faces of a TI which minimizes cross-effects. 
The TI configuration allows us to derive the "mode chan- 
nels" just by the application of the "pattern matrix" 
B. This matrix contains the numerical value of the five 
quadrupolar spherical harmonics at each transducer po- 
sition. So the pattern matrix projects the displacement 
of each transducer to the normal mode components. In 
this way the amplitude of the normal quadrupolar modes 
(the mode channels) can be experimentally monitored by 
the set of transducers. These mode channel amplitudes 
are important quantities in data analysis once they are 
directly related to the tensorial components of the GW 
as we will show later. 

The TI configuration implies Equation to present 17 
second order differential equations (5 from the sphere and 
12 from the set of transducers) which can be solved under 
their symmetry properties. In the TI configuration the 
last six components of the vector w(cc;) correspond to the 
displacement of the last transducer resonator q^^ (^^^ (^f^_ 
nal mass). We can rewrite these last vector components 
in the Equation [3 as 

q^Hc.)=PHFH, (3) 

where P{u;) is the last six lines from the matrix product 
in Equation 121 When we multiply B by Equation [31 we 
obtain an expression for the mode channels which is given 

by 

The quantities denoted by F'^{uj) correspond to driven 
forces which are acting on: the sphere (6*), the inter- 
mediate mass {Ri) and the final mass {R2) The 
transfer function matrix ^{cj) is the proportionality ra- 
tio between the mode channels and the mode effective 
forces. We can easily refer the intermediate and final 
mass noises into the sphere by using the mode response 
function ft^^{uj) to the mode i of the transducers. In this 
way, we projected all noise sources to the antenna modes 
which is measured as 

g{u,) = i{io)F-'\io) (5) 

where F^^\ij) denotes signal plus noise forces which are 
acting on the modes. 

For the moment we can suppose that only a GW sig- 
nal is exciting the detector and see how this interaction 
works. 
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FIG. 2: The Euler angles which define the relative orientation 
between xyz and XY Z. Rotate xyz counterclockwise around 
its axis by a to give xyz. Rotate xyz counterclockwise 
around its y' axis by f5 to give x"y" z" . And finally rotate 
x"y" z" counterclockwise around its z" axis by 7 to give XY Z. 



Such an assumption (a = 0) makes the inverse prob- 
lem easier to solve considering that only two orientation 
angles are important: (5 and 7. These angles easily de- 
termine the incoming direction of a candidate signal if it 
is admitted that x-axis matches with local north-south 
direction and z-axis coincides with local zenith. In this 
case f5 is the complementary angle for the local elevation 
d and 7 is opposite to the azimuth angle 0. If these an- 
gles and the exact arrival time of a signal are known then 
the sky position of the source can easily be found. 

So we can suppose an incident gravitational plane wave 
(which is carrying the reference frame XY which in- 
teracts with the detector (in the reference frame xyz) 
such GW components (in tranverse and traceless form) 
is represented by 







hx{t) -h. 




(9) 



Such components of strain matrix H'(t) would be pro- 
jected into a new strain matrix H(t) by defining H(t) = 
RH'(t)R'''. The projected strain matrix H(t) presents 
components as 17J 



B. The Effective Gravitational Force and the 
Spherical Amplitudes 

Once a passing gravitational wave generates an acceler- 
ation field such effective force (in the frequency domain) 
is given by 



(6) 



we can deconvolve the mode channels in the Equation 
into the spherical amplitudes 



V 



V3' 

h2(t) 



h2{t) hi{t) 
-hi{t) - ^h^{t) hsit) 



(10) 

These components are linear combinations of the five 
components of the spherical amplitudes h presented in 
Equation [3 which are denoted by 

hi (t) = /i+ (t) ^ (1 + cos^ P) cos 27 + /ix {t) cos p sin 27, 

(11a) 



h(/) = 



(7) 



These spherical amplitudes depend on each of the po- 
larization amplitudes and the relative orientation to the 
lab coordinate system (usually centered in the sphere's 
center of mass) [13, [l3 • They can be understood as the 
projected wave polarizations into the sphere's quadrupo- 
lar modes under incident angles a, /3 e 7. 

The relative orientation between two orthogonal 
righthanded 3D cartesian coordinate systems xyz (de- 
tector frame) and XYZ (wave frame) is described by a 
3x3 rotation matrix R, which is commonly parameter- 
ized by three so-called Euler angles a, e 7. Figure [21 
shows the rotation between these two reference frames. 

We can set the detector frame of such form so that 

y = y' 



: y" by setting a = 0, which causes the rotation 



matrix to be expressed by 

cos 7 cos sin 7 

R = I — sin 7 cos jS cos 7 sin 7 sin (3 

sin [3 cos (3 



(8) 



h2{t) 



-h+{t)-{l + cos^ (3) sin 27- 



hx{t) cos /3 cos 27, 
(lib) 



h3{t) = —h-^{t)- sin 2/3 sin7 + hx (t) sin/3cos7, (11c) 



/i4(t) = hj^{t) - sin 2/3 cos 7 -\- hx{t) sin /3 sin 7, (Hd) 



h^it) = h^{t)^V3 sin^ p. (lie) 

These equations show us that the spherical amplitudes 
are extremely dependent on GW polarizations and inci- 
dent angles. With these assumptions one observes that 
they also can mimetize both signal and noise depending 
on what kind of source has been considered. 
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Sciienber]g Sensitivity at +.2 K 
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FIG. 3: The SCHENBERG sensitivity curves at 4.2 K (upper 
line) for a single quadrupole mode and the individual contri- 
butions of the noise sources. 



C. The Noise Sources and the Sensitivity Curve 

Actually as one can see the effective gravitational force 
in Equations Q to [3 is not the only force that acts on 
the sphere. The resultant force ¥^^^{00) has other com- 
ponents that come from noise sources which are acting 
on the whole system (e.g. thermal Langevin's forces). 
Taking into account that we are focused on generating a 
data time series similar to that which comes from a real 
detector these sources must be taken into consideration 
when simulating because they will determine the detector 
behavior. 

Many solutions to minimize the noise sources contri- 
butions have been studied and a lot of what has been 
learned over the last decades with bar instruments has 
been reused and reorganized to be applied to spherical 
detectors (e.g. vibration isolation systems, cryogenic de- 
vices, electromagnetic isolation, etc). But some noise 
sources are intrinsic to the instrument and hard to elim- 
inate. 

We are interested in the detector sensitivity and its 
optimum bandwidth when usual noises (thermal, back- 
action, series and phase noises) [l^ are acting on it. In 
order to determine how large the sensitive bandwidth 
Afm is for the mode m we have adopted the relation 



/ + CX) 
-CO 



\hm (^) I'^duj, 



(12) 



where /i^-^ represents the minimum squared value of 
hni{(jo). The knowledge of spectral densities of each of 
the noise sources makes it possible to determine their 
contributions to the mode channels and consequently the 
values of hm {(jo). 

Figure 01 presents the sensitivity curves for a 
quadrupole mode (all others are considered the same. 



which they actually are for the degenerated case) and 
the individual contributions from inputted noise sources. 
These noise sources are referred into the sphere via trans- 
fer function (Equations from ^ to 0). 

This fact permits the location of the GW source in 
the sky, as well as the determination of the polarization 
amplitudes of the GW. 



D. The Signal plus Noise 

To simulate the SCHENBERG behavior when it is ex- 
cited both by an incident gravitational wave and noise we 
needed to use some waveform. We adopted the Astro- 
GravS Catalog ( |http : //a str ogra vs . gsf c . nasa . govj ) 
where a large number of gravitational waveforms are 
available. As an example, in this paper we used an 
adiabatic inspiral of a corotating binary which was 
rescaled to double 1.9 M© neutron stars (a chirp signal). 

Both signal and noise are expressed in distinct units 
(i.e. at the frequency domain the signal amplitude is ex- 
pressed in Hz~^ and the noise amplitude is expressed in 
, SO we turn signal amplitude into another impor- 
tant quantity the characteristic strain which is defined as 



hc{f) = y A/|/i^^(/)P, where \h^^ {f)\^ is the power 
spectral density of the signal at frequency / El^. The 
characteristic strain is essentially the rms signal in a fre- 
quency interval of width A/ centered at frequency / and 
it has same units of noise sources. If the noise floor is 
flat then this assumption can be used without restrictions 
otherwise a correction must be applied. Our correction 
function has a shape that looks like the "inverse" of the 
sensitivity curve. It assines the right weight for each 
point of the sensitivity curve. So, the points of maxi- 
mum sensitivity are multiplied by one and, for example, 
a point in the curve with sensitivity 100 times worse is 
multiplied by 1/100. 

After this correction we can introduce signal plus noise 
to model and estimate the signal-to-noise ratio needed for 
detection. 



E. The Amplitude Signal-to-Noise Ratio (ASNR) 

The amplitude signal-to-noise ratio (ASNR) was esti- 
mated from Thorne (1987) 22] as 



ASNR 



W{f)^-^^^^^df (13) 



where W{f) represents the applied filter function. SN{f) 
is the equivalent noise profile. It denotes the sum of the 
noises which acts on each normal mode. It is given by 



(14) 



m=l 
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where hm{f) is the sensitivity of mode m as presented in 
Figure [31 It is easy to show that this sum represents the 
total noise energy on the sphere. This relation takes in 
account all sky sensitivity because 

f2hl^h' = hl + h%. (15) 

m=l 

In other words it has no incident direction dependence. 
As we presumed at the beginning the spherical detector 
really is omnidirectional. 
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The Simulated Data 
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In order to generate time series which contain the am- 
plitudes w(t) and, consequently, q2 (the displacement of 
the second mechanical mode of the transducers), we in- 
troduced first only known noise sources. 

We admitted that noise forces behave like the sum of 
sinusoids with random phases at frequencies — oo < a; < 
oo. However it is impossible to sum an infinite num- 
ber of sinusoids. The Sampling Theorem guarantees that 
only frequencies less than the sampling frequency (we use 
/sample = 16384 Hz) could be detected. So we simulated 
at the frequency domain the sinusoids with amplitudes 
Ai{ij) = ^J Si{(jo) that corresponds to spectral densities 
of noise forces at frequency uo and we gave them random 
phases so that 



F{uo) — Ai (lo) (cos (p -\- i sin cp) , 



(16) 



where <(/:?< 27r is the chosen phase (with flat dis- 
tribution) and i = is the complex unit. The In- 
verse Fourier Transform sums the sinusoids to generate 
the time series. Acting in this way we obtained the noise 
forces at the time domain. These values are inputted into 
the model. 

The data present a gaussian distribution as expected. 
With this data to hand we can apply the inverse problem 
to obtain the spherical amplitudes and find an incident 
position for each sample. Figure ^ presents a histogram 
of the spherical amplitudes which are obtained from the 
set of six transducers. The out of band contribution was 
eliminated by applying a rectangular band-pass filter (the 
width of such filter is obtained by Equation E| and it 
must be centered in the resonant frequency fm for the m 
mode). 

From Figure^one can also observe that sometimes one 
mode is more excited {hs in this collection of data) than 
another, but the degree of excitation changes from one 
moment to another which reflects the system dynamics. 
The energy passes from one mode to another via trans- 
ducers and also because broken symmetry of the sphere 
due to the hole machined for the suspension system. 

With these quantities we are able to find the "incident 
direction" and the "GW polarization amplitude" for each 
sample by the inverse problem solution. 



FIG. 4: Histogram of event values for the five spherical am- 
plitudes hm- Their distribuition is perfectly gaussian for a 
noise-only case. 



III. THE INVERSE PROBLEM SOLUTION 

In order to obtain H'(t) from H(t) only the angles /3 
and 7 are needed once we adopt the a = convention. 
Equation EH can be multiplied by from left and by R 
from right and one obtains 

R^H(t)R = R^RH'(t)R^R. (17a) 

Once RtR = RRt = I we have 

H'(t) = R^H(t)R. (17b) 

The numeric values of /3 and 7 can be easily obtained by 
imposing that the last column of H'(t) is null (it really 
is according to GR, GW has no effect on propagation 
direction z'). Therefore the solutions has the following 
format 



tan 7 
tan/3 



= ± 



y 

1 

X 

y 1 

z sin 7 



(18) 
(19) 



The minus sign in Equation ^| comes from the Euler y- 
convention and the double sign in Equation spots the 
fact of the sphere is incapable of distinguishing antipodal 
sources. 

The H'(t) determinant is null and the equation system 
is indeterminate - one of the equations depends on the 
others. However only two of them are needed for the uni- 
vocal determination of 13 and 7. At this time we impose 
a contour condition such that when tan 7 = we have 
tan/3 = 00, just as it should be 23]. In this way we have 



2^3 



(20a) 



2^3 



v/3 



hi{t)h^{t) + ^h^{tY - hs{ty, (20b) 
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h2{t)hs{t) + h^{t) hi{t) + V^sW . (20c) 



we take into account the H'(t) aspect {hj^ = = —H'22 
and = = ^21)- One can observe that H(t) and 
H'(t) are symmetrical so the values of ^21 
always identical and noise independent. However this 
restriction is not applicable to and H22 so we take 
In order to obtain an expression for GW polarization = (i^n + ^22)/2 |l3- So, in this way, we found that 

I 




hj^{t) = (l + cos^/3) ( -hi{t)cos2j- -h2{t)sm2j +sin2/3 — hsit) sin 7 + - (t) cos 7 + — /^sWsin^/?, (21a) 



hx{t) = COS P {hi (t) sin 27 + /12 (^) cos 27) + sin p {hs (t) cos 7 + /14 (^) sin 7) . (21b) 

I 



These equations can also be derived from Equation [TTl 
and represent the solution of the inverse problem. They 
are valid either for noisy and noiseless cases. The noise 
distribution over the sphere surface can now be estimated 
for each sample once "a direction" is obtained. 

IV. THE PROBLEM OF INCOMPLETENESS 

The fact that only the quadrupolar modes have been 
considered instead of all spherical harmonic functions 
leads to a problem of incompleteness in the sphere's re- 
construction. 




FIG. 5: The sum of the five quadrupolar modes over the 
sphere in a Hammer- Aitoff projection (the brightest areas rep- 
resent the highest absolute amplitudes). 

The sum of the five quadrupolar modes (which is shown 
in Figure with equal degrees of excitation implies in 
amplified regions which points to a preferential direction 
(Figure when only noise is considered. This direction 
is actually perpendicular to the position of higher ampli- 
tudes of sphere vibrations. 

This non-uniform distribution leads us to produce a 




FIG. 6: When the direction is calculated the pattern points 
to some preferential direction. So, the noise isn't uniformly 
distributed over the sphere. 



"gain mask" which normalizes the output position map. 
We call it "gain mask" because either signal or noise that 
incides in that direction is more easily detected by the 
inverse solution once it excites all sphere modes equally. 
An analytic expression for this is not easily obtained so 
we use a numerical method by using a large number of 
noise only simulations (~ 100 of 32 s at 16384 Hz) to make 
it smooth enough to apply it to the problem. 

The gain mask enhances areas of low noise levels in 
the output position map in order to produce an uniform 
distribution of noise for all points in the map. We intro- 
duced a signal such as the one mentioned in Subsection 
III Dl and located the source at different distances to dis- 
cover what ASNR is needed to detect it. 



V. THE RESULTS 

Once the original source position (in the simulated 
data) is known we can estimate the precision on its de- 
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termination. A source must appear at position (^0,^0) 
like a point if no noise is affecting the detector, however 
in the case of noise it could be found at another place. 
Low signal-to-noise ratios cause an indetermination in 
the source location and results in an error. This error 
puts the source in a place in a region of the sky Q (a 
solid angle) delimited by AO and Acj) around its original 
position. Such an error can be quantified by the solid an- 
gle ft ^ Joq-ao 1^0- OdOdcj). The larger the solid 
angle filled by the source the larger the error on the po- 
sition determination is. This effect is presented in Figure 



observer 




FIG. 7: The solid angle described by a source if we have a 
noisy detector. The found position could be inside a solid 
angle delimited by and Acj) around its original position 
(^0, 0o)- Such an area determines the error in position recov- 
ering. 




FIG. 8: Position of a source with ASNR ^ 3. This figure 
represents the mean of ^ 150 distinct simulations. 

Figure IHl shows the incident direction of a source with 
ASNR ~ 3 after we had used a simple band-pass filter 
on the data. The figure represents the mean amplitude 
in ~ 150 distinct simulations. It is easy to notice that a 
brighter area appears after these number of simulations. 



Such an area describes a solid angle that corresponds to 
the error on the source location. All pointed directions 
by inverse solution are inside this area. So the area of 
the source location for ASNR ~ 3 case is equivalent 
Vt ^ 0.235 sr which approximately corresponds to 2% of 
the whole sky. 




FIG. 9: Position of a source with ASNR ^ 3. This figure 
represents the signal in only one simulation. 

However, in a single simulation, as shown in Figure 
ini it is already possible to detect the source position. 
It fills not just a point in the sky but something like a 
strip during the time that the chirp signal remains in the 
sensitivity bandwidth. On the other hand, noise samples 
mark only points randomly distributed over the sphere 
surface. So behavior like the one showed in Figure [HI is a 
candidate signature. 

6e-22 I 1 1 1 1 1 1 1 



4e-22 




-6e-22 I 1 1 1 1 1 1 1 

5000 10000 15000 20000 25000 30000 

sample number 

FIG. 10: Pseudo-amplitudes from the inverse solution. The 
dots represents 3cr of the sample distribution. 

This ASNR ~ 3 represents a trustable limit for a can- 
didate detection. As one can see in Figure 11171 the level of 
the signal is almost diving in the noise but its maximum 
amplitude is > 3cr which makes the detection possible. 
The search for signals with lower values for ASNR needs 
another approach, a statistical one. 
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VI. CONCLUSIONS 

We numerically solved the inverse problem to recover 
the information about an input GW source such as its 
incident direction and polarizations. It is possible to de- 
termine event energy and to find a direction for its source 
even for ASNR ~ 3. The amazing result emerges from 
the fact that a spherical detector emulates 5 indepen- 
dent detectors (the mode channels) with an important 
peculiarity: the simultaneous correlation between them. 
Neither interpolation nor another artifice is needed. 

The mode channels come directly from a product be- 
tween the pattern matrix and the data and such a com- 
putational operation is easily done in real time. How- 
ever, the spherical amplitudes (Equation [TJ depends on 
the inverse of the transfer function matrix ^{co) which 
is difficult to calculate. A solution comes from the fact 
that such quantity does not depend either on signal or on 
noise. It only depends on the electro-mechanical features 
of the detector, so it could be stored on disk or even on 
the computer memory to be used in a real time analy- 
sis. So the presented method can be applied to a real 
time data analysis software which will be able to emit an 
on-line warning to be checked. 

This method is very fiexible and can be adapted for the 
real case. Although the model is highly parametrized so 



the values of the physical parameters must be known. 
However a set of transfer function with a distinct col- 
lection of parameters could be tested offline to validate 
some candidates. 

With this method it becomes evident that a spherical 
GW detector, like the Mario SCHENBERG, wih be able 
to detect nearby chirp events (up to Magellanic clouds) 
at its first goal sensitivity (4.2 K) just by using a simple 
band-pass filter. 

For ASNR < 3 (using a simple band-pass filter) an- 
other approach must be used. We tested the efficiency 
of a matched filter and checked that it works well when 
the waveform and incident direction are known. However 
it has a poor performance when the input signal doesn't 
match perfectly with the template. But this will be in- 
vestigated in detail in a future work. 
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